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Abstract 

The  orientation  of  patches  on  the  surface  of  an  object  can  be  determined  from  multiple 
images  taken  with  different  illumination,  but  from  the  same  viewing  position  This  method, 
referred  to  as  photometric  stereo,  can  be  implemented  using  table  lookup  based  on 
numerical  inversion  of  experimentally  determined  reflectance  maps.  Here  we  concentrate  on 
objects  with  specularly  reflecting  surfaces,  since  these  are  of  importance  in  industrial 
applications.  Previous  methods,  intended  for  diffusely  reflecting  surfaces,  employed  point 
source  illumination,  which  is  quite  unsuitable  in  this  case.  Instead,  we  use  a  distributed 
light  source  obtained  by  uneven  illumination  of  a  diffusely  reflecting  planar  surface. 
Experimental  results  are  shown  to  verify  analytic  expressions  obtained  for  a  method 
employing  three  light  source  distributions.  - 


This  report  describes  research  done  at  the  Artificial  Intelligence  Laboratory  of  the 
Massachusetts  Institute  of  Technology.  Support  for  the  Laboratory's  artificial  intelligence 
research  is  provided  in  part  by  the  Office  of  Naval  Research  under  contract  N000H-77-C- 
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Notation 

©  :  zenith  angle  (local  coordinate  system) 

:  azimuth  angle  (local  coordinate  system) 

6  :  zenith  angle  (viewer  coordinate  system) 

<f>  :  azimuth  angle  (viewer  coordinate  system) 
F  :  flux  [  watt/sr  ] 

E  :  irradiance  [  watt/m2  ] 

L  :  radiance  (  watt/(m2sr)  ] 

--  suffix  -- 
g  :  source 
i :  incident 
e  :  emitting 
p  :  image  (or  picture) 

„ :  surface  normal 
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0  Overview 


The  technique  of  photometric  stereo  lias  been  proven  by  an  experiment  whose  schema  is 
shown  in  Fig.  I  Orientations  of  surface  patches  on  a  sphere  obtained  using  this  technique 
are  shown  in  Fig.  2. 


kax-search 


Fig.  I  The  overall  schema  of  the  experiment.  The  experiment  has  two  parts:  initialization, 
represented  by  broken  lines,  and  the  lookup  of  surface  orientations  for  a  given  object  using 
measured  image  brightness.  Brightness  arrays  obtained  from  a  TV  camera  are  normalized. 
A  lookup  table  is  constructed  in  the  initial  process  using  reflectance  map  techniques  and 
Newton's  method.  The  table  is  used  to  obtain  surface  orientations. 
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Fig.  2  Needle  diagram  showing  (projected)  surface  normals  on  the  surface  of  a  metallic 
sphere.  Each  needle  shown  sticking  up  out  of  a  point  represents  the  orientation  of  a  surface 
patch 
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I  Introduction 

The  geometric  d  endence  of  the  reflectance  characteristics  of  a  surface  can  be  expressed  in 
terms  of  the  direction  from  the  surface  to  the  viewer,  the  direction  from  the  surface  to  the 
light-source  and  the  surface  normal  The  apparent  brightness  (scene  radiance)  seen  by  the 
viewer  can  be  expressed  as  a  function  of  the  three  angles  between  these  directions 

Alternatively,  we  can  express  the  dependence  in  terms  of  the  slope  components  p  and  q, 
used  as  axes  in  gradient  space  [I],  There  are  several  ways  to  project  the  unit  surface  normal 
onto  a  plane  In  the  "traditional"  method  one  projects  a  point  on  the  unit  sphere  onto  a 
tangent  plane  from  a  center  placed  at  one  pole  of  the  unit  sphere  The  plane  is  tangent  to 
the  sphere  at  the  opposite  pole  (stereo-graphic  projection).  This  projection  may  be  useful 
when  we  consider  characteristics  of  a  surface  which  has  both  specularity  and  transparency 
like  the  surface  between  air  and  water  or  windows  of  the  John  Hancock  building;  one 
hemisphere  corresponds  to  reflected  brightness  distribution,  the  other  to  transmitted 
brightness  distribution. 

p’»  0z/0x/(  yu(0z/dx)z«(8z/0y)z  t  I  )  (  I  ) 

q'=  0z/0y/(  ,/l*(0z/0x)2*(0z/0y)2  t  I), 

where  z=Z(x,y)  and  the  sign  is  chosen  depending  on  whether  the  point  is  in  the  upward 
hemisphere  or  the  downward  hemisphere  In  another  method,  which  has  been  used  by 
Horn  [I],  one  projects  the  unit  surface  normal  onto  the  tangent  plane  from  the  cento  of  the 
unit  sphere  (central  projection). 

p  =  Az/Ax  (  2  ) 

q  -  Az/Ay 


Fig  ’•  I  loin's  projection  and  the  stereographic  projection  onto  a  plane 


Although  the  tiaditional  method  may  be  more  familiar  and  allows  one  to  express  any 
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dnection  without  ambiguity,  the  zenith  angle  is  projected  into  half  of  the  zenith  angle  and. 
thus,  the  coordinate  transformation  is  not  as  convenient  as  it  is  for  the  second  method. 
Moreover,  we  only  need  to  consider  the  upper  hemisphere  of  possible  directions,  since 
surfaces  turned  away  from  the  viewer  cannot  be  seen  in  any  case.  It  is  for  these  reasons 
that  we  have  chosen  the  gradient  space  projection. 

If  we  take  the  direction  from  the  surface  to  the  viewer  as  the  direction  of  the  z-axis,  then 
the  reflectance  properties  of  a  surface  patch  depend  on  (p,q),  the  direction  of  the  surface 
normal  and  (ps.qs),  the  direction  to  the  source  [I],  Each  point  in  gradient  space  (the  tangent 
plane),  corresponds  to  a  particular  surface  orientation  (based  on  the  direction  of  the  viewer). 
If  we  know  the  reflectance  characteristics  of  an  object,  we  can  calculate  how  bright  a  surface 
element  with  that  orientation  will  appear.  It  is  convenient  to  use  contour  lines  to  connect 
those  points  in  gradient  space  which  correspond  to  surface  orientations  which  give  rise  to 
the  same  apparent  brightness.  It  is  because  of  these  contour  lines  that  the  resulting  diagram 
is  referred  to  as  the  "reflectance  map"  [I],  It  is  denoted  by  R(  p,q  ). 

Using  the  reflectance  map,  the  basic  imaging  equation  is 

E,(x,y)  -  R,(  p.q.x.y  ),  (  3  ) 

where  E,(x,y)  is  the  brightness  (image  irradiance)  in  the  image-forming  system  at  the  point 
(x.y)  in  the  image  plane  This  equation  contains  two  unknown  variables  p,  q  and  one 
quantity  E,.  which  can  be  measured  in  the  image. 

In  the  above  equation,  the  subscript  j  is  used  to  denote  different  illumination  conditions, 
for  each  value  of  the  subscript,  a  different  image  is  obtained,  and  a  different  reflectance 
map  applies  If  two  images  are  taken,  two  such  equations  provide  constraints  on  the 
possible  values  of  p  and  q  This  permits  us  to  solve  for  the  gradient.  Because  of  the  non¬ 
linearity  of  the  equations,  however,  a  number  of  solutions  may  be  found  at  times.  In  this 
case,  a  third  image  may  be  used  to  disambiguate  the  remaining  possibilities.  This  is  the 
pi  inciple  of  photometric  stereo  [2.3,1 

Ortliogiaphic  projection  can  simplify  the  calculation  considerably.  If  we  can  assume  that 
the  object  is  small  compared  with  the  distance  to  the  source  and  the  image  forming  system, 
then  the  viewer  direction  can  be  approximated  as  the  axis  of  the  image-forming  system,  and 
we  can  ti rat  t he  system  as  orthographic.  There  are  two  merits  to  this  approximation.  One 
is  that  we  can  neglect  the  effect  of  position.  The  right  hand  side  of  Eq.  3  depends  only  on 
(p.q),  namely,  we  can  apply  the  same  reflectance  map  on  all  point  in  the  image  Another 
benefit  is  that  we  can  calculate  R(  p,q  )  more  easily  because  the  approximation  fixes  the 
viewing  direction.  So  we  can,  for  example,  rotate  the  source  keeping  the  phase  angle 
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constant  (The  phase  angle  is  the  angle  between  the  source  and  the  viewer,  measured  at  the 
object).  This  means  that  we  can  obtain  a  new  reflectance  map  just  by  rotating  the  old  one 
[?]. 

2  Rasic  Consideration 

2.1  Relationship  between  Radiance  and  Irradiance 

One  of  the  main  points  of  our  discussion  here  is  that  we  consider  only  specular  components 
of  reflectance  when  we  calculate  the  reflectance  map,  since  many  industrial  materials  are 
marie  of  metal  and  have  strong  specularity  and  little  diffuse  reflection  Experiments  show 
that  only  I  or  27.  of  the  incident  light  is  reflected  diffusely  from  some  metallic  surfaces,  with 
most  of  the  rest  reflected  specularly  We  cannot  treat  this  kind  of  material  using  the  usual 
Lambertian  model  for  reflection  of  light  from  a  surface.  It  is  also  clearly  inappropriate  to 
use  point  sources  to  illuminate  such  a  surface,  since  very  few  surface  patches  will  be  oriented 
coirectly  to  reflect  any  light  and  we  will  only  see  virtual  images  of  the  point  sources 

We  consider  two  kinds  of  light  sources.  One  is  an  extended  light  source.  The  other  is  a 
collimated  one.  An  example  of  an  extended  light  source  is  a  fluorescent  light  fixture.  A 
collimated  light  may  be  generated  by  using  a  laser  and  an  inverted  telescope. 

Three  relationships  exist  between  a  light  source  and  the  image  plane  The  first  one  is  that 
between  source  radiance  and  incident  irradiance  on  a  sample  surface.  Next,  and  most 
important,  is  that  between  incident  irradiance  and  emitting  radiance  from  the  surface  (this 
is  captured  in  the  reflectance  map).  The  last  one  is  that  between  emitting  radiance  and 
image  irradiance.  Since  the  two  kinds  of  light  sources  have  somewhat  different 
relationships,  we  consider  them  separately. 

(extended  source) 

The  relationship  between  source  radiance  from  an  extended  light  and  the  incident 
irradiance  is  fairly  straightforward.  We  assume  a  surface  patch  of  a  light  source  (cl A)  is 
emitting  energy  to  a  surface  patch  of  a  sample  surface  (dB).  The  emitted  energy  from  the 
source  incident  on  the  sample  surface  is 

dFs  =  Lsdw,dAcosO|  (  *  ) 

■=  Ls(l/rz)dAdBcos©8cos0t. 

where  <.>,  is  a  solid  angle  corresponding  to  the  surface  patch  dB  viewed  from  dA,  and  r  is 
the  distance  between  the  source  and  the  sample  surface.  On  the  other  hand,  the  received 
flux  is 
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Fig.  4  Relationships  between  light  source,  surface,  and  image  plane  are  depicted. 


dF,  =  L,dwsd  BcosO,  (5) 

»  L,(l/r2)dAdBcos©iCos@g, 

where  ws  is  a  solid  angle  corresponding  to  the  surface  dA  of  the  source  viewed  from  dB 

If  L.s  and  L,  can  be  regarded  as  constant  within  the  small  areas  dA,  dB,  then 

L,  -  W  (  6  ) 

This  implies  that  we  can  regard  the  "incident  radiance"  as  the  source  radiance  in  the  same 
direction. 

For  a  specular  surface  and  an  extended  source,  we  can  obtain  from  the  definition  of  the  bi¬ 
directional  reflectance  distribution  function  (BRDF)  [1,4] 

L.(©..+.)-//f.,idEi 

■=  //26(sin29,-sin2©,)  L,(©r*,)cos0^in©ld©j  d'F, 


-  L,(0t, *,♦»), 

(7) 

wheie 

fp ,  --  2<Xsinz©t-sinz©,)6(4,#-4,,*7r) 

(8) 

is  the  BRDF  for  a  perfectly  specular  surface;  L|(9f.*,or)  is  the  "incident  radiance",  and 
is  the  reflected  radiance.  This  result  is  also  predicted  from  elementary  physical 
laws.  Thus,  even  though  the  source  distribution  may  be  complicated,  only  the  contribution 
from  a  single  direction  (©„♦,»*)  need  be  considered  and  one  need  not  be  concerned  with 
the  effects  of  other  areas  of  the  distributed  light  source  Strictly  speaking,  L*  is  not  defined 
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as  a  function  of  local  angles  but  viewer  angles  Light  source  which  we  can  construct  has  a 
fixed  brightness  distribution  on  the  viewer  system.  However,  it  is  easier  to  understand  the 
result  expressed  using  the  above  notation.  In  the  next  section  we  will  give  the  equations  for 
the  transformation  between  the  viewer  coordinate  system  and  the  local  coordinate  system 

The  relationship  between  reflected  (scene)  radiance  and  image  irradiance  is 

Lp  *-  { ( 7r  /  4  )(d  /fp)2cos4« }  L,,  (  9  ) 

where  f„,  d,  «  are  effective  focal  length  of  the  lens,  the  diameter  of  the  entrance  apeiture 
and  the  off  axis  angle,  respectively!  2  ]. 

Finally,  using  on  Eq.  6,7,9  we  see  that  image  irradiance  at  a  particular  point  Is  proportional 
to  the  source  radiance  in  a  direction  which  depends  on  the  orientation  of  the  corresponding 
surface  patch!  I  ].  That  is,  the  brightness  of  a  particular  surface  patch  is  simply  equal  to  the 
brightness  of  the  part  of  the  extended  source  which  it  happens  to  reflect.  Thus,  even 
though  the  source  distribution  may  be  complicated,  only  the  contribution  from  a  single 
direction  (pf,qs)  need  be  considered  at  any  time. 

(collimated  light  source) 

In  the  case  of  a  collimated  source,  it  is  easy  to  use  irradiance  to  measure  relationship 
between  source  intensity  and  incident  brightness  Clearly,  provided  that  a  surface  patch  is 
small  compared  to  the  diameter  of  the  source,  source  irradiance  Es  is  the  same  as  E, 
measured  on  a  surface  lying  orthogonal  to  the  rays  According  to  [  I  ],  incident  radiance  L, 
is 

L,(©„4',')  =  E,6(©,-©s)6('f',-+s)/sin©s.  (  |Q  ) 

The  beauty  of  this  expression  is  that  the  double-delta  expression  has  dimension 
[  l/steradian  ].  because  when  it  is  integrated  with  respect  to  solid  angle  the  result  has  no 
dimension.  .So  the  total  dimension  of  the  right  hand  side  of  Eq.  10  is  [  watt/(m2sr)  ]  and  it 
coi responds  to  the  left  hand  side.  If  all  the  direction  are  covered  with  collimated  sources 
each  of  which  has  a  particular  brightness  distinguished  with  the  direction, 

L, <©„*,)  =  EkE.^Q,  -©Sk)6(4>l-4>sk)/sin0n  (  II  ) 

where  k  denotes  the  k-th  light  source. 

Now  we  can  calculate  the  scene  radiance  by  using  the  BRDF[  1,4  ].  The  BRDF  for  the 
sin  fare  is  the  same  as  in  the  former  case  and  is  given  in  Eq  8 

■  // Lj  (0„4’,)26(sin20,-sin20j)6(4t-4>j*ir)cos0isin0ld0Id4,j  (  |2  ) 

-  1,(0,,  4yir) 

Again,  the  same  abbreviation  is  used;  namely,  brightness  distribution  Ek  is  not  defined  as 
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a  function  of  the  local  coordinates  This  will  be  discussed  in  the  next  section. 

The  relationship  of  scene  radiance  and  image  irradiance  is,  of  course,  the  same  as  for  the 
extended  light  source.  In  this  case  we  also  find  that,  using  Eqs  II,  12,  8,  the  image 
irradiance  can  be  calculated  from  the  source  irradiance  determined  by  the  orientation  of  the 
local  normal  on  the  surface  patch. 

(transformation  from  a  local  coordinate  system  to  a  viewer  coordinate  system) 

It  is  convenient  to  change  from  the  local  coordinate  system  to  the  viewer-centered  coordinate 
system  [1]  because  the  reflectance  map  is  defined  in  a  viewer-centered  coordinate  system,  and 
the  light  source  distribution  is  given  based  on  the  direction  of  the  viewer. 

In  both  cases,  we  have  to  determine  the  viewer  coordinates  corresponding  to 
when  the  viewer  and  the  surface  normal  have  the  relationship  specified  by  (©e,4e).  By 
using  the  relations 
©,  -  ©e  =  en, 

4’i  =  4’,,  *  ir,  (  13  ) 

we  get  the  corresponding  directions  of  the  viewer  coordinate 

e,  =  2  en 

<t>,  =  </>„.  (  H  ) 

This  relationship  is  obtained  either  by  considering  Fig.  5  or  substituting  Eq  13  into 
spherical  trigonometric  relations. 

Thus,  we  can  rewrite  Eq.  7  and  Eq.  12  using  the  viewer-center  coordinate  system. 

Ep(0n.<f>„)  =  «L(0n,</>„)  *  aL(20n, <*>„), 
where  ct  is  a  constant.  We  observe  this  Ep  as  image  irradiance 

Fiom  Eq.  14  we  can  also  obtain 
Ps  -  2  p„  /  (l  pnz-qn2) 

qs  -  2  q„  /  (l-Pn  V>  <  >5 »  ) 

We  can  finally  express  brightness  distribution  in  the  gradient  space.  In  the  case  of  an 
extended  light  source,  there  is  no  problem. 

R<Pn-qn)  ■  L.(p8q8).  (  16.1  ) 

wheie 

pn  -  dz/dx  =  -cos<f>ntan0„ 
q„  =  dild y  -  -sin</>ntan0„. 


Fig  5  Relationships  between  a  local  coordinate  system  and  a  viewer  coordinate  system 


If  we  use  a  collimated  light  source, 

L,(20n.<f\,>  =  ZkFk/26{0n-esk/2)^n-(|)9k)/sin0sK 

=  Es(2#n,</>„)/(2sin  20n), 

where  Es  is  obtained  from  k  such  that  Psk  -2  6n  and  <f>,k  »  4>n.  We  finally  get 
R(pn,qn)=  Es(p9.qs)/(4(l.p8zcq,2)2). 


r 
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2.2  Consideration  on  Types  of  Brightness  Distribution. 

The  mapping  (pr,.qn)eU-{(x,y)|(x2*y2)<l}  to  (pt>q,)eV-R 2-{oo}  has  the  following 
characteristics 

(i>  (pn-q„)=o  ->  (ps.q$M) 

(2)  pn-ft  ->  ps-0 

(3)  q„=0  ->  qs=Q 

<*0  Pn/qn'Ps/qs  (qn.qs  are  non-zero) 

(5)  p2*qz=(yi*psz*q,z  -  0Z/(ps2*qs2) 

Actually, the  forward  transformation  is  given  in  Eq.  15.1.  Namely, 

Ps  -  2  p„  /  0-pZ-q2) 

qs  =  2  qn  /  (l-p2-q2).  (  15 1  ) 

On  the  other  hand,  the  inverse  transformation  is 

Pn » ps</hS7-  0/(ps2*qs2) 

qn  -  qs</i*ps2tq«2  -  D/(p,2*qs2)  ( <5  2  ) 

We  will  use  a  differentiable  single-valued  function  as  the  reflectance  function  Roughly 
speaking,  since  we  determine  surface  orientation  from  brightness,  If  the  function  is  not 
single  valued,  the  inverse  function  theorem  does  not  apply,  and  this  makes  the  situation 
very  difficult 

From  the  consideration  given  in  the  Appendix,  the  reflectance  map  Is  a  regular  surface 
based  on  the  choice  of  the  reflectance  function  as  a  differentiable  single-valued  function  It 
means  that  at  each  point  on  the  surface  each  point  possesses  a  neighborhood  where  we  can 
define  an  infinite  derivative  on  the  surface  independent  of  choices  of  parameterization  It 
aho  turns  out  that  (p„,qn)  and  (pt,qs)  are  two  parameterizations  of  the  surface  This  means 
that  each  point  on  the  surface  uniquely  corresponds  to  a  particular  point  (p„.qft)  and  (ps.qs) 
and  that  (pn,qn)  and  (pt,qs)  are  representations  of  the  point  by  using  two  different  coordinate 
systems  Wc  can  treat  the  surface  using  one  or  the  other  coordinate  at  our  convenience  and 
we  can  always  tianslate  from  one  to  the  other  without  ambiguity.  These  properties  will  be 
used  when  wc  build  the  lookup  table  using  Newton's  method. 

2  3  Shape  of  Light  Source  and  its  Brightness  Distribution 

Basically,  we  considered  two  shape  for  making  a  light  source  under  the  condition  of 
differentiable  one-to-one  mapping.  One  is  a  sphere  and  the  other  Is  a  plane.  Table  I 
contains  their  shapes  and  brightness  distributions. 


lkeuchi  fc  Horn 


Photometric  Stereo 


Ikeuchi  &  Horn 


13 


Photometric  Stereo 


KAX-SHArC!! 


Fig  6  The  overall  schema  of  the  experiment.  The  technique  consists  of  two  types  of  jobs: 
off-line  preparation  of  tables  (represented  by  broken  lines)  and  on-line  processing.  Image 
brightness  is  obtained  by  using  a  TV  camera.  We  search  for  the  maximum  value  of 
brightness  in  each  array.  Brightness  arrays  are  calibrated  and  normalized.  Surface 
orientations  are  obtained  from  the  lookup  table. 

3.2  Consideration  of  Light  Source 

We  use  a  Lambertian  surface  as  a  source  plane.  It  is  illuminated  by  a  linear  lamp  as  shown 
in  Fig  7.  Though  a  spherical  shaped  source  can  easily  cover  directions  of  more  than  ninety 
degrees,  it  is  difficult  to  build  such  a  device  and  difficult  to  control  the  distribution  of  the 
light  on  it.  particularly  if  internection  are  taken  into  account.  On  the  other  hand,  if  we  use 
a  planar  source  illuminated  by  a  lamp,  the  brightness  distribution  is  complicated,  but  we 
can  calculate  a  reflectance  map  as  in  case  2  of  Table  I.  It  is  possible  to  cover  angles  of  more 
than  ninety  degrees  by  making  a  box-like  source.  In  that  case,  however,  we  have  to  treat 
each  plane  separately  because  the  surface  normal  is  not  differentiable  at  the  intersection  of 
two  planes.  Thus,  we  consider  one  plane  surface  which  is  assumed  to  have  the  Lambertian 
characteristics  in  order  to  obtain  some  analytic  results  to  help  in  the  design  of  the  system. 
This  plane  is  illuminated  by  a  line  source  as  a  representative  case. 

Brightness  distribution  on  the  surface  is  calculated  using  Eq.  18.  Let  f  be  the  flux  late  per 
unit  souice  length  (  watt/(m  sr)].  Then  the  total  irradiance  E  [  watt/mz  J  is 

E  -  /fcos0|Cosfl2/rzdt.  (  17  ) 

From  cos#,  -  yxz»l0z/r  and  cos02  -  Iq/t.  we  can  finally  get 
E  -  (f10/2az)[ltan'‘(y«L)/a  -  tan''(y-L)/a) 

♦|  a(y*L)/(az*(y»L)z)-a(y-L)/(az*(y-L)z) }  J,  (  18  ) 
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Fig.  7  The  object  is  illuminated  solely  by  the  light  reflected  of  an  overhead  Lambertian 
surface.  This  surface  receives  light  from  a  linear  lamp  positioned  below  so  as  not  to 
illuminate  the  object.  The  TV  camera  peers  through  a  small  hole  in  the  overhead  surface. 
Brightness  distribution  from  a  linear  lamp  on  the  Lambertian  surface  can  be  calculated  in 
Eq  18  We  used  three  linear  lamps,  placed  symmetrically,  120  degree.  In  the  center  of  the 
lamps  there  is  a  hole  The  object  is  observed  through  the  hole. 


Ikeuchi  &  Horn 


15 


Photometric  Stereo 


where  a  =  J*1* l0^,  L  is  the  length  of  the  line  source,  and  1q  is  the  distance  between  the  line 
source  and  the  Lambertian  surface. 

We  put  an  object  just  under  the  point  (x,y)  «  (0.5I(>,0.0).  E(x,y)  is  a  symmetric  function  with 
respect  to  y.  Thus,  it  is  natural  to  put  an  object  somewhere  on  the  x  axis.  E(x,0)  is  an  s- 
shape  function  and  has  a  inflection  point  near  x«0.5lo.  At  that  point 
E(0f>lo,0,0)=07l5E(0. 0,0.0).  So  it  is  convenient  to  put  an  object  just  under  this  point.  If  we 
denote  the  distance  from  the  surface  to  the  object  as  l(, 
x  -  •lPstW2 

y  -  *iq*  <  w ) 


Finally,  we  can  get  a  reflectance  map 

R<pft.qn)  -  [flo/2{<-llP.*lo/2)2ViX 

tan'l{(-l,qs*lo/2V^-l,pt*l(>/2)z*l0z)-tan- '  {(-I  ,q$-lo/2)/ 

♦i  V^ip8*lo/2)z*l0z(-llqt*l0/2)}/{(-llpt*l0/2)z*l0z*(-l|q<«l()/2)zi 

-{ >X‘iIp«410/2)z«loZ('llq*'l0/2))/{{-llpt*l0/2)z*l02*(-l|qs-lo/2)2)l  (  20  ) 

where  p,  and  q,  are  as  defined  above,  and  L  -  Iq,2I|  -  Iq  because  these  values  are  found  to 
be  optimal  upon  simulation. 

We  used  three  linear  lamps,  placed  symmetrically,  120^  apart,  about  the  hole  through  which 
the  object  is  observed.  One  lamp  is  turned  on  at  a  time,  giving  rise  to  a  reflectance  map 
like  the  one  shown  in  Eq.  20. 

By  using  three  reflectance  maps,  we  can  determine  the  surface  orientation.  Theoretically  it 
is  possible  to  determine  the  orientation  by  using  two  maps  only.  However,  there  may  be 
more  than  one  solution  because  of  the  non-linearity  of  the  equations.  If  we  use  three  maps, 
we  can  determine  a  unique  solution  easily  as  shown  in  Fig.  8. 

3.3  llow  to  Make  a  Lookup  Table 

The  most  convenient  method  for  converting  a  triple  of  measured  brightness  to  an 
orientation  is  by  means  of  a  lookup  table  made  from  the  reflectance  maps.  This  table, 
indexed  by  quantized  brightness  measurements,  contains  surface  orientation.  Although  it  is 
possible  to  make  a  three  dimensional  lookup  table  in  which  each  dimension  corresponds  to 
the  brightness  of  a  surface  patch  under  one  of  the  three  sources,  the  weakest  brightness 
contains  relatively  large  measurement  errors  and  we  use  it  only  to  choose  a  solution  between 
two  alternatives.  Thus,  the  lookup  table  can  be  two  dimensional.  We  look  up  an  entry 
using  the  two  largest  brightness  values.  Each  entry  contains  two  alternative  solutions.  Each 


CELL  INTENSITY  . 

RTJ&BER  SOURECE  1! '  2  3 

183  25?  e. $17849997  0.232034993  0.404375553 

108  253  0.53948421  0.477330793  0.83733455  A  [ 

103  259  0.69183O9S  0.496798813  0.852332614 

103260  e. 78222797  .0.58294154  0.91585983  |b  1 

108  261  0.81491274  !  0.59471041  0.94813397  |  j  P=  -0.1 

103  262  0.75035409  I  0.63671693  0.89895891  J  Q=_0.4 

103268  0.O4396324&  0.8522275  0-701 1291  |L 

108  269  8.78222797  10.83536255  0.518473156  * 

108  2?e  0.51946413  10.560397836  0.264884242  D 

Fig  8  How  to  use  the  reflectance  map  method  Three  reflectance  maps  drawn  in  the  lame 
gradient  space.  Since  the  lamps  are  symmetrically  configured,  brightness  distribution 
coriesponding  to  a  linear  lamp  is  also  symmetric  having  the  same  shape  I  hus,  wc  need  to 
calculate  only  one  distribution  and  then  rotate  it  120  degrees.  The  resulting  function  is  the 
desired  one 

We  can  determine  (p.q)  from  three  values  of  brightness,  where  each  brightness  value 
tonesponds  to  each  source  condition  For  example,  switching  on  souice  1  will  yield  a 
btightness  value  of  0  7822  at  cell  B  Similarly,  source  2  and  source  3  yield  0S829  and  0.91  VI 
respectively  hem  the  above  diagram,  we  note  that  the  point  (  -0.1  -04  )  satisfies  this  triple. 
Hence,  the  surface  orientation  at  the  cell  B  is  (  01  -0.4  ). 
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alternative  solution  also  contains  a  surface  orientation  and  the  third  brightness 

We  construct  the  lookup  table  by  using  Newton’s  method.  An  elegant  method  exists  in  the 
case  of  a  Lambertian  reflector  and  point  sources  [21  In  our  case,  however,  we  have  to 
calculate  it  numerically  Specifically,  we  solve  two  expressions  like  Eq  20  for  p  and  q  The 
expressions  differ  in  the  brightness  distribution  of  the  source,  obtained  by  turning  on  one  of 
the  three  linear  lamps 

There  are  two  choices  for  the  calculation.  One  method,  which  we  rejected,  divides  the 
gradient  space  into  a  mesh  and  gets  a  solution  by  interpolation  A  similar  technique  is 
considered  in  [  6  ]  to  make  an  exact  color  from  a  combination  of  three  primal  colors  T  hey 
tried  to  find  the  nearest  grid  point  to  a  solution  and  then  determine  which  neighboring  box 
(the  small  space  between  grid  points)  contains  the  solution.  Finally,  they  used  an 
inteipolation  method  from  vertices  of  the  box  to  determine  the  result.  The  technique  would 
be  useful  when  a  reflectance  map  is  obtained  experimentally. 

The  other  method  is  to  solve  Eq  20  numerically  for  p  and  q  In  case  that  reflectance  maps 
ate  given  in  analytic  form,  it  makes  sense  to  use  Newton’s  method.  Moreover,  we  have 
found  that  Newton’s  method  converges  within  four  iterations  with  relative  error  smaller 
than  0  001,  provided  that  the  initial  point  is  well  chosen.  Since  the  neighboring  solutions  of 
the  lookup  table  are  good  initial  points,  the  calculation  requires  the  order  of  4N2,  where  N 
is  the  number  of  mesh  points,  while  the  former  method  needs  at  least  N4/2  calculations 
besides  the  initial  calculations  at  each  mesh. 

The  N-dirnensional  Newton  method  involves  solving  the  simultaneous  linear  equations: 

.  i"dfI/dxklxll  »  x^'W")  -  0.  (  21  ) 

where  f,  is  the  i-th  element  of  a  vector  function  of  n  dimensions  and  xj1'  is  the  i-th  element 
of  a  solution  vector  of  the  j-th  iteration. 

From  the  consideration  of  Section  2.2,  we  can  treat  the  reflectance  map  as  a  function  of 
(p,.qs)  when  we  get  a  solution  using  Newton’s  method  Newton’s  method  finds  an 
intei  sections  of  two  surfaces  Since  Eq  20  is  differentiable,  the  result  is  a  regular  surface 
and  the  itnation  conveiges  to  a  point  aiong  the  regular  surface  From  the  consideration  in 
Appendix  and  Section  2  2,  this  situation  can  be  observed  from  either  of  two 
parameter izations  and  the  expression  of  the  result  in  one  parameterization  should 
corie«pond  to  the  other  exactly.  Thus,  we  use  a  function  parameterized  by  (pt.qt)  ns  f,  we 
obtain  a  numerical  solution  and  then  we  project  the  solution  to  (p„.q„)  space. 
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We  can  rewrite  Eq  21  as  Eq.  23.  in  our  case  of  a  planar  source  illuminated  by  a  linear  lamp. 
The  two  mapping  functions  are 

f<ivqs) "  R(ps.qs)'Ep,n 

g(ps.qs)  =  RfptCosa-qjSina.pjSina^qsCOsaVEpn,,  (  22  ) 

where  a  is  an  angle  between  two  light  sources  and  Ein,  E,m  are  image  brightness 
coi  responding  to  the  (n,m)  element  of  the  lookup  table  We  can  finally  get 

xf"‘U  x^-fFV'F.  <  23  ) 

where  X  =  '(p,.qs),  F  *  ’(f.g)  and  (F )‘*  is  a  inverse  of  the  Jaccobian  matrix 

The  complete  algorithm  is  shown  in  Fig  9.  At  first,  we  try  to  find  solution  corresponding  to 
a  brightness  pair  <0  4.0  4)  from  initial  points  (10,1.0)  and  (-1.0, 1.0).  Next,  we  proceed  to  get 
solutions  along  a  line  on  which  the  first  brightness  is  constant  (0.4),  using  former  solutions 
as  new  initial  points  until  two  corresponding  solutions  from  two  different  directions  are  the 
same  oi  we  cannot  get  solutions  within  five  iterations.  Then,  we  go  to  the  next  mesh  point 
along  a  line  on  which  the  second  brightness  is  constant  (0.4).  After  the  entire  area  is 
processed,  the  algorithm  terminates. 


Fig  9  How  to  get  entire  solutions  using  Newton’s  method. 
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3.3  Mow  to  Get  a  Surface  Orientation 

Image  bnghtness  is  obtained  from  the  TV  camera.  To  reduce  the  noise  typical  of  these 
devices,  we  took  more  than  one  picture  per  light  source,  and  arrays  corresponding  to  the 
same  light  source  were  averaged.  The  resulting  three  brightness  arrays,  one  for  each  light 
source,  were  the  input  to  the  photometric  stereo  system. 

We  search  each  array  for  its  maximum  brightness.  If  objects  are  convex  and  the  visual 
angle  is  wide  enough,  the  image  array  always  contains  a  brightness  value  corresponding  to 
the  maximum  source  brightness.  By  normalizing  the  results  using  these  maximum  values, 
we  can  cancel  the  effect  of  varying  albedo  and  image  sensor  response.  Namely,  we  can 
identify  the  maximum  value  in  an  image  array  with  the  value  1.0  in  the  reflectance  map. 
However,  we  cannot  normalize  the  values  of  the  image  arrays  at  this  stage  because  of  the 
non  linearity  of  the  TV  camera  Normalization  will  be  done  after  the  value  is  calibrated. 
In  other  word,  since  the  output  of  the  TV  camera  is  not  linearly  proportional  to  a 
brightness  value,  we  have  to  normalize  the  output  after  we  convert  the  output  to  a  certain 
value  proportional  to  image  brightness.  It  is  interesting  that  the  maximum  point  always 
contains  important  information  and  this  search  resembles  Land’s  lightness  search  model 
[71 

Biightness  calibration  is  done  using  a  six  step  KODAK  gray  scale.  Normalization  is  based 
on  the  maximum  value  So  each  point  on  the  image  has  a  triple  of  calibrated  brightness 
under  thiee  different  light  sources. 

From  the  lookup  table  we  can  get  the  surface  orientation.  The  lookup  table  entry  is 
accessed  using  the  two  largest  brightness  values.  The  lookup  table  is  actually  two  triangular 
two  dimensional  matrices.  To  increase  the  accuracy  of  computation  and  economy  of  memory 
use,  the  first  dimension  represents  the  largest  image  brightness.  The  second  dimension 
represents  to  the  second  largest.  Each  element  of  the  matrix  contains  the  coriesponding 
suiface  orientation  and  the  smallest  image  brightness. 

Since  the  non-linearity  gives  rise  to  two  solutions  for  each  intensity  pair,  we  have  to  decide 
winch  one  is  the  real  one.  We  choose  a  surface  orientation  by  comparing  the  distance 
between  the  actual  third  image  brightness  and  the  element  of  the  matrix.  The  third 
brightness  is  always  weak  and  likely  to  contain  errors.  It  is  only  used  to  decide  which  of  the 
two  alternatives  obtained  from  the  other  brightness  Is  the  correct  one 
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/  4  Experiment  and  Discussion 

t 

Oui  experimental  results  are  shown  in  Fig.  10.  Fig.  10.1.1  shows  three  brightness  arrays 
The<e  arrays  are  input  information  to  the  photometric  system.  No  suiface  normals  ate 
shown  in  aieas  where  insufficient  information  was  available  in  the  three  images  to 
dcretmme  them  accurately  (see  Fig.  101.2).  The  choice  of  light  source  distribution  a  fleets  the 
extent  of  these  regions  as  well  as  the  accuracy  with  which  surface  normals  can  be  found 


Fig  10  1  I  Three  brightness  arrays.  Each  array  corresponds  to  one  of  the  three  lamps. 
These  arrays  are  input  information  to  the  photometric  stereo  system. 
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Fig  101.2  Direct  output  from  the  system  No  surface  normals  are  shown  in  areas  where 

thei p  was  insufficient  information 


Ikpuchi  &  Horn 


21 


Photometric  Steieo 


We  meet  two  kinds  of  relaxation  methods  The  first  one  is  suitable  to  on  line  systems  This 
method  reads  from  the  lookup  table  twice  exchanging  the  second  and  the  third  brightness 
values  and  determines  the  solution  by  averaging  the  two  previous  solutions  This  simple 
scheme  can  give  good  results  as  shown  in  Fig.  10.1.3.  The  other  method  is  to  find  solutions 
which  minimize  the  difference  of  the  actual  brightness  and  the  theoretical  brightness  with  a 
constraint  derived  from  surface  continuity  [8]  Theoretical  brightness  is  calculated  ficm 
onentation  and  reflectance  maps  Actually,  this  method  is  done  iteratively  The  assumption 
of  sui  face  continuity  requires  that  neighboring  points  should  project  to  neighboring  points 
on  the  gaussian  sphere.  Fig  1014  shows  the  result  obtained  using  the  iterative  method. 
You  can  notice  that  over  a  very  wide  area  the  method  can  find  surface  orientations  This 
method  is  appropriate  for  off-line  systems  Fig  10 1.5  is  a  generated  surface  fiom  their 
sin  face  orientations 
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Fig  10.1.3  Output  from  the  first  relaxation 
method  This  method  would  be  suitable  for  a 
hand -eye  system  where  real-time  response  is 
vital.  The  system  reads  from  the  lookup  table 
twice,  exchanging  the  second  and  the  third 
brightness  values  The  solution  is  determined 
by  averaging  the  two  previous  values.  This 
simple  schema  can  give  good  results. 
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Fig  10 1  4  Output  from  the  iterative  relaxation 
method  This  method  finds  solutions  which 
minimize  the  difference  of  the  actual 
brightness  and  the  theoretical  brightness 
using  a  constraint  derived  from  surface 
continuity. 


(b)  A  surface  along  the  x  axis.  The  bold  line 
represents  the  actual  surface  and  the  diagram 
is  obtained  from  the  photometric  stereo 


(c)  A  surface  along  the  y-axis. 


Fig.  10.1.5  Generated  surfaces  from  the  needle  diagram. 
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A  direct  application  of  our  technique  is  an  industrial  hand-eye  system  that  picks  up  an 
object  out  of  a  jumble  of  material  (see  Fig  10.2).  Although  our  technique  cannot  correctly 
determine  the  surface  orientation  when  there  is  mutual  illumination,  it  does  provide  the 
means  for  detecting  this  condition,  since  the  three  measurements  will  be  inconsistent  In  this 
fashion  erroneous  results  are  avoided.  This  is  important,  since  the  manipulator  might 
otherwise  be  sent  to  a  position  where  it  would  collide  with  other  parts. 
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Fig.  10.2.1  Binary  image  of  a  jumble  of  nuts 
and  bolts. 


Fig.  10.2.2  Needle  diagram  obtained  from 
image  of  nuts  and  bolts. 


Fig  10.2.3  The  generated  surfaces  are  elevated  in  order  to  easily  distinguish  from  the 
surrounding  area. 
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Another  application  of  this  technique  is  the  inspection  of  the  surface  condition  of  metals  If 
a  surface  has  a  crack,  stain,  or  finger  print,  the  image  brightness  triple  yields  inconsistent 
values  in  the  area  of  the  blemish  (see  Fig.  10.3). 


Fig  10  3.1  A  picture  of  a  hook  which  contains  a  lot  of  cracks  and  stains. 
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Fig  10.32  Needle  diagram  of  the  hook  There  exists  a  printed  figure  ”3"  in  the  central  area 
of  the  object.  This  can  be  seen  as  randomness  of  needles  in  the  area.  At  the  low  right  side 
the  surface  is  rather  dark  and  the  system  is  unable  to  determine  the  surface  orientations  at 
that  area  Note  that  a  relaxation  method  would  fail  to  identify  this  central  region  as 
different  from  its  surroundings 
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This  technique  may  be  combined  with  the  Marr-Poggio-Grimson  stereo  technique  now 
being  developed  [  9,10  ].  Their  technique  detects  depth  cues  and  works  well  when  the  image 
contains  many  discontinuities.  On  the  other  hand  our  method  detects  surface  oncntation 
directly  and  works  well  when  the  object  is  smooth.  The  combination  can  be  established  if 
the  output  of  one  of  the  two  stereo  cameras  is  fed  to  our  system  while  the  two  outputs  from 
two  camera  ape  fed  to  their  system.  We  feel  that  the  composite  system  will  produce  an 
excellent  representation  of  the  object  just  as  a  people  are  believed  to  use  both  stereo  and 
shading  information  to  construct  a  symbolic  image  of  the  visual  world. 

Acknowledgment 

The  authors  would  like  to  extend  their  sincere  appreciation  to  Prof.  P.  H.  Winston  of  MIT. 
Discussions  with  H.  Murota  of  the  University  of  Tokyo  and  R.  Sjoberg,  W.  Silver  of  MIT 
were  very  helpful.  Thanks  go  to  B.  Roberts,  K.  Stevens  of  MIT  and  M.  Brooks  of 
University  of  Essex  for  proofreadings. 


References 

[I]  B.  K.  P.  Horn  and  R.  W.  Sjoberg:  "Calculating  the  Reflectance  Map,"  Applied 
Optics.  Vol  18,  No.  11,1979. 

12]  R  J.  Woodham:  "Photometric  Stereo:  A  reflectance  map  technique  for  determining 
surface  orientation  from  Image  Intensity,"  Proc.  S.P.I.E  .Vol.  155, 1978. 

[3]  B  K.  P  Horn,  R.  J.  Woodham,  and  W.  Silver:  "  Determining  Shape  and  Reflectance 
Using  Multiple  Images,"  Al-Memo  490, 1978,  AI  Lab,  MIT. 

[4]  F.  E.  Nicodemus  et  al:  "  Geometrical  Considerations  and  Nomenclature  for 
Reflectance,  "  NBS  Monograph  160, 1977,  U  S.  Department  of  Commerce. 

[5]  M.  P.  D.  Carmo:  "  Differential  Geometry  of  Curves  and  Surfaces."  1976,  Prentice- 
Hall,  New  Jersey. 

[6]  D.  M  Erway:  "  Exact  Color  Reproduction,"  B.S.  thesis,  1978,  EECS.  MIT. 

[7]  E.  Land  and  J  McCann:  "  Lightness  and  Retinex  Theory,"  J  Optical  Society.  Vol.  61, 
No.  1. 1971. 

[8]  T.  Strat:  "  A  numerical  Method  for  Shape  from  Shading  from  a  Single  image,"  M  S. 
thesis.  1979.  EECS,  MIT. 

[9]  D.  Marr  and  T.  Poggio:  "  Cooperative  Computation  of  Stereo  Disparity"  AI-Memo 
364, 1976,  Al  Lab,  MIT. 

[10]  W  E.  L.  Crimson  and  D.  Marr:  "  A  Computer  Implementation  of  a  Theory  of 
Human  Stereo  Vision,"  (to  be  published) 


Ikruchi  &  Horn 


26 


Photometric  Stereo 


Appendix 

Tim  appendix  has  two  purposes;  one  is  to  discover  what  kind  of  characteristics,  grneially 
speaking:,  aie  appropriate  for  a  reflectance  function,  and  the  other  is  to  answer  to  these 
questions  about  the  solution  of  non-linear  equations  in  one  space  is  same  as  that  which  is 
obtained  by  solving  corresponding  equations  in  another  space  and  being  transfoimed  with 
a  conformal  mapping. 

Proposition  I 

A  surface  (ps,qs.  R(ps.qs))  is  a  regular  surface  based  on  the  choice  of  R(p8,q8)  such  that  RZ  U 
->  R  is  a  diffeientiable  one-value  function. 


Fig.  A.i 


(pi  oof) 

A  tegular  surface  satisfies 

I  Its  mapping  function  X  is  differentiable 

2.  X  is  a  homeomorphism 

3.  X  satisfies  the  regularity  condition!  i  1 

whet  p  we  can  denote  X(p8,q8Mp8,q8,  R(p8,q8))  as  the  mapping  function  from  a  open  set  U  in 
R2  to  a  open  set  W  in  R3 

It  is  obvious  that  X  is  differentiable  with  respect  to  p8  and  q8,  because  of  the 
differentiability  of  R(ps.qs).  Condition  3  requires  at  least  one  of  the  three  two-dimensional 
minot  jaccobtan  is  always  non-zero.  Since  cKp8,q$)/d(p$.q8)  is  always  one,  condition  3  is 
satisfied  To  each  unique  (p8.q8)  there  exists  R(p8.q8)  which  is  a  one-to-one  corresponding  to 
(p„.q8)  because  R(ps,qs)  is  one-value  differential  function.  This  implies  X  is  one-to-one 
diffetential  mapping.  From  satisfaction  of  condition  3  and  this  one-to-one  mapping,  inverse 
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theoiem  guarantees  the  existence  of  differential  inverse  X'1.  So  X  is  a 
hompomorphism.  (qed) 

Now  we  ha  vp  known  that  the  surface  is  a  regular  surface  provided  that  we  choose  a  single¬ 
valued  differentiable  function  as  a  reflectance  function;  the  surface  has  various 
characteristics  as  not  that  dependent  on  the  choice  of  the  parameterization  but  intrinsic 
ones  The  following  proposition  is  rather  obvious. 

Proposition  2  The  reflectance  map  of  a  specular  surface  can  be  parameterized  with  (pn.q„) 
such  that  Eq  A  I,  provided  that  we  choose  a  single-valued  differentiable  function  as  a 
reflectance  function 

P»  =  2  P „  /  0-Pn2V) 

qs  =  2  q„  /  (l-p„Z-q„2)  <  A  1  > 


(proof) 

A  mapping  Y:VcRz->UcR2  such  that  (p,.q,)  -  Yfp^qjis  a  differentiable  map,  where  U.  V 
are  open  sets  in  R2.  From  the  theorem  of  the  chain  rule  of  maps,  Z(pB,qnMpn.qn.R(Pnqn)) 
aho  a  differentiable  map.  (qed) 

For  the  leaders  convenience: 

Inverse  Function  tlieorem[5] 

Let  F:UcR°->Rft  be  a  differentiable  mapping  and  suppose  that  at  PcU  the  diffeiential 
dFp  R"->Rn  is  an  isomorphism  Then,  there  exists  a  neighborhood  V  of  p  in  U  and  a 
neighborhood  W  of  F(p)  in  R"  such  that  F:v->W  has  a  differentiable  inverse  F*1  W->V 
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The  Chain  Rule  for  Maps[5] 

Let  F  UcRn->Rm  and  G:VcR"’->R1'  be  differentiable  maps,  where  U  and  V  are  open  sets 
such  that  F(U)cV.  Then,  G°F:U->Rk  is  a  differentiable  map,  and 
d<G*F)p=dGF(p)#dFp.  peU. 


